clear;
close all;
paths;

load('soln.mat');
load('data.mat');
data.yr1_hs= infl_e(2,:)'; data.yr1_cd= infl_e(3,:)';

% unemployment gap
ue = readtable("unemp.csv"); ue.YYYYQ = year(ue.DATE)*10 + quarter(ue.DATE);
ue.uegap = ue.UNRATE - ue.NROUST;
ue = join(data, ue);

px1 = readtable("agg_px.csv");  weight = readtable("weight.csv"); 
ue.yr1_avg = (ue.yr1_hs .* weight.WT1 + ue.yr1_cd .* weight.WT3)./(weight.WT1 + weight.WT3);

% oil price
oilprice = readtable("OILPRICE.csv"); oilprice.YYYYQ = year(oilprice.DATE)*10 + quarter(oilprice.DATE);

%% correlation with oilprice

oilprice = join(oilprice, data);
soln.corr_cd = corr(oilprice.px1_educg3 - oilprice.yr1_cd, oilprice.OILPRICE);
soln.corr_hs = corr(oilprice.px1_educg1 - oilprice.yr1_hs, oilprice.OILPRICE);

%% slope of Phillips Curve

recession = 20074;

pc_sample1 = ue.YYYYQ >= 19851 & ue.YYYYQ < recession;
pc_sample2 = ue.YYYYQ >= recession & ue.YYYYQ < 20222;
pc_sample3 = ue.YYYYQ >=19791;


x1 = linspace(min(ue.uegap), max(ue.uegap), 1000);
c1_1 = polyfit(ue.uegap(pc_sample1), ue.inflation(pc_sample1) - ue.px1_educg1(pc_sample1),1);
c1_2 = polyfit(ue.uegap(pc_sample2), ue.inflation(pc_sample2) - ue.px1_educg1(pc_sample2),1);
z1_1 = polyval(c1_1, x1); z1_2 = polyval(c1_2, x1);

pc1_nc = fitlm(ue.uegap(pc_sample1), ue.inflation(pc_sample1) - ue.px1_educg1(pc_sample1), 'Intercept', true);
pc2_nc = fitlm(ue.uegap(pc_sample2), ue.inflation(pc_sample2) - ue.px1_educg1(pc_sample2), 'Intercept', true);

figure;
scatter(ue.uegap(pc_sample1), ue.inflation(pc_sample1) - ue.px1_educg1(pc_sample1), '*', 'MarkerEdgeColor',[0.3010 0.7450 0.9330], 'LineWidth', 1); hold on; 
plot(x1, z1_1, '-', 'Color', [0.3010 0.7450 0.9330], 'LineWidth', 2); 
scatter(ue.uegap(pc_sample2), ue.inflation(pc_sample2) - ue.px1_educg1(pc_sample2), 'o', 'MarkerEdgeColor',[0 0.4470 0.7410], 'LineWidth', 1);
plot(x1, z1_2, '--', 'Color', [0 0.4470 0.7410], 'LineWidth', 2);
legend('1985Q1-2007Q3',' ','2007Q4-2022Q2',' ');
xlabel('unemployment gap','fontsize', 14);
ylabel('\pi - \pi^e: raw non-college', 'fontsize', 16);fontsize(gcf,scale=1.5);
ylim([-6,3])
style = hgexport('factorystyle'); style.Color = 'gray';
hgexport(gcf,"figures/pc_raw_nc.eps",style);


x1 = linspace(min(ue.uegap), max(ue.uegap), 1000);
c1_1 = polyfit(ue.uegap(pc_sample1), ue.inflation(pc_sample1) - ue.px1_educg3(pc_sample1),1);
c1_2 = polyfit(ue.uegap(pc_sample2), ue.inflation(pc_sample2) - ue.px1_educg3(pc_sample2),1);
z1_1 = polyval(c1_1, x1); z1_2 = polyval(c1_2, x1);

pc1_cd = fitlm(ue.uegap(pc_sample1), ue.inflation(pc_sample1) - ue.px1_educg3(pc_sample1), 'Intercept', true);
pc2_cd = fitlm(ue.uegap(pc_sample2), ue.inflation(pc_sample2) - ue.px1_educg3(pc_sample2), 'Intercept', true);

figure;
scatter(ue.uegap(pc_sample1), ue.inflation(pc_sample1) - ue.px1_educg3(pc_sample1), '*', 'MarkerEdgeColor',[0.9290 0.6940 0.1250], 'LineWidth', 1); hold on; 
plot(x1, z1_1, '-', 'Color', [0.9290 0.6940 0.1250], 'LineWidth', 2); 
scatter(ue.uegap(pc_sample2), ue.inflation(pc_sample2) - ue.px1_educg3(pc_sample2), 'o', 'MarkerEdgeColor',[0.8500 0.3250 0.0980], 'LineWidth', 1);
plot(x1, z1_2, '--', 'Color', [0.8500 0.3250 0.0980], 'LineWidth', 2);
legend('1985Q1-2007Q3',' ','2007Q4-2022Q2',' ');
xlabel('unemployment gap','fontsize', 14);
ylabel('\pi - \pi^e: raw college', 'fontsize', 16);fontsize(gcf,scale=1.5);
ylim([-6,3])
style = hgexport('factorystyle'); style.Color = 'gray';
hgexport(gcf,"figures/pc_raw_cd.eps",style);


x1 = linspace(min(ue.uegap), max(ue.uegap), 1000);
c1_1 = polyfit(ue.uegap(pc_sample1), ue.inflation(pc_sample1) - ue.yr1_hs(pc_sample1),1);
c1_2 = polyfit(ue.uegap(pc_sample2), ue.inflation(pc_sample2) - ue.yr1_hs(pc_sample2),1);
z1_1 = polyval(c1_1, x1); z1_2 = polyval(c1_2, x1);

pc1_nc_m = fitlm(ue.uegap(pc_sample1), ue.inflation(pc_sample1) - ue.yr1_hs(pc_sample1), 'Intercept', true);
pc2_nc_m = fitlm(ue.uegap(pc_sample2), ue.inflation(pc_sample2) - ue.yr1_hs(pc_sample2), 'Intercept', true);

figure;
scatter(ue.uegap(pc_sample1), ue.inflation(pc_sample1) - ue.yr1_hs(pc_sample1), '*', 'MarkerEdgeColor',[0.3010 0.7450 0.9330], 'LineWidth', 1); hold on; 
plot(x1, z1_1, '-', 'Color', [0.3010 0.7450 0.9330], 'LineWidth', 2); 
scatter(ue.uegap(pc_sample2), ue.inflation(pc_sample2) - ue.yr1_hs(pc_sample2), 'o', 'MarkerEdgeColor',[0 0.4470 0.7410], 'LineWidth', 1);
plot(x1, z1_2, '--', 'Color', [0 0.4470 0.7410], 'LineWidth', 2);
legend('1985Q1-2007Q3',' ','2007Q4-2022Q2',' ');
xlabel('unemployment gap','fontsize', 14);
ylabel('\pi - \pi^e: purified non-college', 'fontsize', 16);fontsize(gcf,scale=1.5);
ylim([-6,3])
style = hgexport('factorystyle'); style.Color = 'gray';
hgexport(gcf,"figures/pc_pure_nc.eps",style);


x1 = linspace(min(ue.uegap), max(ue.uegap), 1000);
c1_1 = polyfit(ue.uegap(pc_sample1), ue.inflation(pc_sample1) - ue.yr1_cd(pc_sample1),1);
c1_2 = polyfit(ue.uegap(pc_sample2), ue.inflation(pc_sample2) - ue.yr1_cd(pc_sample2),1);
z1_1 = polyval(c1_1, x1); z1_2 = polyval(c1_2, x1);

pc1_cd_m = fitlm(ue.uegap(pc_sample1), ue.inflation(pc_sample1) - ue.yr1_cd(pc_sample1), 'Intercept', true);
pc2_cd_m = fitlm(ue.uegap(pc_sample2), ue.inflation(pc_sample2) - ue.yr1_cd(pc_sample2), 'Intercept', true);

figure;
scatter(ue.uegap(pc_sample1), ue.inflation(pc_sample1) - ue.yr1_cd(pc_sample1), '*', 'MarkerEdgeColor', [0.9290 0.6940 0.1250], 'LineWidth', 1); hold on; 
plot(x1, z1_1, '-', 'Color', [0.9290 0.6940 0.1250], 'LineWidth', 2); 
scatter(ue.uegap(pc_sample2), ue.inflation(pc_sample2) - ue.yr1_cd(pc_sample2), 'o', 'MarkerEdgeColor',[0.8500 0.3250 0.0980], 'LineWidth', 1);
plot(x1, z1_2, '--', 'Color', [0.8500 0.3250 0.0980], 'LineWidth', 2);
legend('1985Q1-2007Q3',' ','2007Q4-2022Q2',' ');
xlabel('unemployment gap','fontsize', 14);
ylabel('\pi - \pi^e: purified college', 'fontsize', 16);fontsize(gcf,scale=1.5);
ylim([-6,3])
style = hgexport('factorystyle'); style.Color = 'gray';
hgexport(gcf,"figures/pc_pure_cd.eps",style);
